{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# FSLE experiment in med\n\nExample to build Finite Size Lyapunov Exponents, parameter values must be adapted for your case.\n\nExample use a method similar to `AVISO flse`_\n\n    https://www.aviso.altimetry.fr/en/data/products/value-added-products/\n    fsle-finite-size-lyapunov-exponents/fsle-description.html\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from matplotlib import pyplot as plt\nfrom numba import njit\nfrom numpy import arange, arctan2, empty, isnan, log, ma, meshgrid, ones, pi, zeros\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.data import get_demo_path\nfrom py_eddy_tracker.dataset.grid import GridCollection, RegularGridDataset\n\nstart_logger().setLevel(\"ERROR\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## ADT in med\n:py:meth:`~py_eddy_tracker.dataset.grid.GridCollection.from_netcdf_cube` method is\nmade for data stores in time cube, you could use also\n:py:meth:`~py_eddy_tracker.dataset.grid.GridCollection.from_netcdf_list` method to\nload data-cube from multiple file.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "c = GridCollection.from_netcdf_cube(\n    get_demo_path(\"dt_med_allsat_phy_l4_2005T2.nc\"),\n    \"longitude\",\n    \"latitude\",\n    \"time\",\n    # To create U/V variable\n    heigth=\"adt\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Methods to compute FSLE\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "@njit(cache=True, fastmath=True)\ndef check_p(x, y, flse, theta, m_set, m, dt, dist_init=0.02, dist_max=0.6):\n    \"\"\"\n    Check if distance between eastern or northern particle to center particle is bigger than `dist_max`\n    \"\"\"\n    nb_p = x.shape[0] // 3\n    delta = dist_max**2\n    for i in range(nb_p):\n        i0 = i * 3\n        i_n = i0 + 1\n        i_e = i0 + 2\n        # If particle already set, we skip\n        if m[i0] or m[i_n] or m[i_e]:\n            continue\n        # Distance with north\n        dxn, dyn = x[i0] - x[i_n], y[i0] - y[i_n]\n        dn = dxn**2 + dyn**2\n        # Distance with east\n        dxe, dye = x[i0] - x[i_e], y[i0] - y[i_e]\n        de = dxe**2 + dye**2\n\n        if dn >= delta or de >= delta:\n            s1 = dn + de\n            at1 = 2 * (dxe * dxn + dye * dyn)\n            at2 = de - dn\n            s2 = ((dxn + dye) ** 2 + (dxe - dyn) ** 2) * (\n                (dxn - dye) ** 2 + (dxe + dyn) ** 2\n            )\n            flse[i] = 1 / (2 * dt) * log(1 / (2 * dist_init**2) * (s1 + s2**0.5))\n            theta[i] = arctan2(at1, at2 + s2) * 180 / pi\n            # To know where value are set\n            m_set[i] = False\n            # To stop particle advection\n            m[i0], m[i_n], m[i_e] = True, True, True\n\n\n@njit(cache=True)\ndef build_triplet(x, y, step=0.02):\n    \"\"\"\n    Triplet building for each position we add east and north point with defined step\n    \"\"\"\n    nb_x = x.shape[0]\n    x_ = empty(nb_x * 3, dtype=x.dtype)\n    y_ = empty(nb_x * 3, dtype=y.dtype)\n    for i in range(nb_x):\n        i0 = i * 3\n        i_n, i_e = i0 + 1, i0 + 2\n        x__, y__ = x[i], y[i]\n        x_[i0], y_[i0] = x__, y__\n        x_[i_n], y_[i_n] = x__, y__ + step\n        x_[i_e], y_[i_e] = x__ + step, y__\n    return x_, y_"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Settings\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Step in degrees for ouput\nstep_grid_out = 1 / 25.0\n# Initial separation in degrees\ndist_init = 1 / 50.0\n# Final separation in degrees\ndist_max = 1 / 5.0\n# Time of start\nt0 = 20268\n# Number of time step by days\ntime_step_by_days = 5\n# Maximal time of advection\n# Here we limit because our data cube cover only 3 month\nnb_days = 85\n# Backward or forward\nbackward = True"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Particles\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x0_, y0_ = -5, 30\nlon_p = arange(x0_, x0_ + 43, step_grid_out)\nlat_p = arange(y0_, y0_ + 16, step_grid_out)\ny0, x0 = meshgrid(lat_p, lon_p)\ngrid_shape = x0.shape\nx0, y0 = x0.reshape(-1), y0.reshape(-1)\n# Identify all particle not on land\nm = ~isnan(c[t0].interp(\"adt\", x0, y0))\nx0, y0 = x0[m], y0[m]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## FSLE\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Array to compute fsle\nfsle = zeros(x0.shape[0], dtype=\"f4\")\ntheta = zeros(x0.shape[0], dtype=\"f4\")\nmask = ones(x0.shape[0], dtype=\"f4\")\nx, y = build_triplet(x0, y0, dist_init)\nused = zeros(x.shape[0], dtype=\"bool\")\n\n# advection generator\nkw = dict(t_init=t0, nb_step=1, backward=backward, mask_particule=used, u_name=\"u\", v_name=\"v\")\np = c.advect(x, y, time_step=86400 / time_step_by_days, **kw)\n\n# We check at each step of advection if particle distance is over `dist_max`\nfor i in range(time_step_by_days * nb_days):\n    t, xt, yt = p.__next__()\n    dt = t / 86400.0 - t0\n    check_p(xt, yt, fsle, theta, mask, used, dt, dist_max=dist_max, dist_init=dist_init)\n\n# Get index with original_position\ni = ((x0 - x0_) / step_grid_out).astype(\"i4\")\nj = ((y0 - y0_) / step_grid_out).astype(\"i4\")\nfsle_ = empty(grid_shape, dtype=\"f4\")\ntheta_ = empty(grid_shape, dtype=\"f4\")\nmask_ = ones(grid_shape, dtype=\"bool\")\nfsle_[i, j] = fsle\ntheta_[i, j] = theta\nmask_[i, j] = mask\n# Create a grid object\nfsle_custom = RegularGridDataset.with_array(\n    coordinates=(\"lon\", \"lat\"),\n    datas=dict(\n        fsle=ma.array(fsle_, mask=mask_),\n        theta=ma.array(theta_, mask=mask_),\n        lon=lon_p,\n        lat=lat_p,\n    ),\n    centered=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Display FSLE\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(13, 5), dpi=150)\nax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\nax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\nax.set_aspect(\"equal\")\nax.set_title(\"Finite size lyapunov exponent\", weight=\"bold\")\nkw = dict(cmap=\"viridis_r\", vmin=-20, vmax=0)\nm = fsle_custom.display(ax, 1 / fsle_custom.grid(\"fsle\"), **kw)\nax.grid()\n_ = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.01, 0.9]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Display Theta\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(13, 5), dpi=150)\nax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\nax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\nax.set_aspect(\"equal\")\nax.set_title(\"Theta from finite size lyapunov exponent\", weight=\"bold\")\nkw = dict(cmap=\"Spectral_r\", vmin=-180, vmax=180)\nm = fsle_custom.display(ax, fsle_custom.grid(\"theta\"), **kw)\nax.grid()\n_ = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.01, 0.9]))"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}